!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


module mod_ti

 implicit none
 public :: c_ode, i_rhs_ev, i_add, i_stimes
 private

 type, abstract :: c_ode
 contains
  procedure(i_rhs_ev), deferred, pass(uuu) :: rhs_ev
  procedure(i_add),    deferred, pass(z) :: add
!  procedure(i_add),    deferred, pass    :: add
!  generic :: operator(+) => add
  procedure(i_stimes), deferred, pass(y) :: stimes
!  procedure(i_stimes), deferred, pass(x) :: stimes
!  generic :: operator(*) => stimes
 end type c_ode

 abstract interface
  pure subroutine i_rhs_ev(tnd,uuu)
   import :: c_ode
   implicit none
   class(c_ode), intent(in)  :: uuu
   class(c_ode), intent(out) :: tnd
  end subroutine i_rhs_ev
 end interface

! abstract interface
!  pure function i_add(x,y) result(z)
!   import :: c_ode
!   class(c_ode), intent(in) :: x, y
!   class(c_ode), allocatable :: z
!  end function i_add
! end interface
!
! abstract interface
!  pure function i_stimes(r,x) result(y)
!   import :: c_ode
!   real, intent(in) :: r
!   class(c_ode), intent(in) :: x
!   class(c_ode), allocatable :: y
!  end function i_stimes
! end interface

 abstract interface
  pure subroutine i_add(z,x,y)
   import :: c_ode
   class(c_ode), intent(in) :: x, y
   class(c_ode), intent(out) :: z
  end subroutine i_add
 end interface

 abstract interface
  pure subroutine i_stimes(y,r,x)
   import :: c_ode
   real, intent(in) :: r
   class(c_ode), intent(in) :: x
   class(c_ode), intent(out) :: y
  end subroutine i_stimes
 end interface

end module mod_ti

!--------------------------------------------------------------------

module mod_ee

 use mod_ti, only: c_ode, i_rhs_ev, i_add, i_stimes

 implicit none
 public :: ee_step

contains

 pure subroutine ee_step(ynew,ynow,dt)
  class(c_ode), intent(in)  :: ynow
  real, intent(in) :: dt
  class(c_ode), intent(out) :: ynew

  class(c_ode), allocatable :: tnd

   !allocate(tnd,source=ynow)
   allocate(tnd,mold=ynow)
   !allocate(class(ynow) :: tnd)

   call ynow%rhs_ev(ynew)
   call tnd%stimes(dt,ynew)
   call ynew%add(ynow,tnd)
   !ynew = ynow + tnd%stimes(dt)
   !ynew = tnd%stimes(dt)
   !ynew = ynow%add(tnd%stimes(dt))
   deallocate(tnd)
 end subroutine ee_step

end module mod_ee

!--------------------------------------------------------------------

module mod_pb

 use mod_ti, only: c_ode, i_rhs_ev

 implicit none
 public :: t_ode, rhs

 type, extends(c_ode) :: t_ode
  real :: y(3)
 contains
  procedure, pass(uuu) :: rhs_ev => rhs
  procedure, pass(z) :: add
  procedure, pass(y) :: stimes
 end type t_ode

contains

 pure subroutine rhs(tnd,uuu)
  class(t_ode), intent(in)  :: uuu
  class(c_ode), intent(out) :: tnd

   select type(tnd)
    type is(t_ode)
     tnd%y = -3.0*uuu%y
   end select
 end subroutine rhs

 pure subroutine add(z,x,y)
  class(c_ode), intent(in) :: x, y
  class(t_ode), intent(out) :: z

   select type(x)
    type is(t_ode)
     select type(y)
      type is(t_ode)
       z%y = x%y + y%y
     end select
   end select
 end subroutine add

 pure subroutine stimes(y,r,x)
  real, intent(in) :: r
  class(c_ode), intent(in) :: x
  class(t_ode), intent(out) :: y

   select type(x)
    type is(t_ode)
     y%y = r*x%y
   end select
 end subroutine stimes

end module mod_pb

!--------------------------------------------------------------------

program ode_test

 use mod_ti, only: c_ode, i_rhs_ev
 use mod_pb, only: t_ode, rhs
 use mod_ee, only: ee_step

 implicit none
 integer :: i
 type(t_ode) :: a, b, tmp

 a%y = (/ 1.0 , 2.0 , 3.0 /)
 do i=1,10
   !call ee_step(b,a,0.1,tmp)
   call ee_step(b,a,0.1)
   write(*,*) b%y
   a = b
 enddo

end program ode_test
